;+
; NAME:
;       OMNI_MASS_FUNCTION
;
; PURPOSE:
;       Computes the mass function based on parameters in the
;       mass-deriv configuration file.
;
; CATEGORY:
;       distance-omnibus Mass Routine
;
; CALLING SEQUENCE:
;       OMNI_MASS_FUNCTION [,CONFFILE=cfile][,NREAL=nreal][,/ML]
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       CONFFILE -- Name of the configuration file to use for survey
;                   information [Default: conffiles/survey_info.conf]
;       NREAL    -- Number of realizations of the mass function to
;                   compute [Default: 1000].
;       ITERN    -- Iteration number for this value of NREAL for the
;                   multiple-realization version.
;       MMAX     -- Maximum mass (in M_sun) to be imposed on the mass
;                   function fitting -- object masses above this value
;                   will not be considered in the fit.  [Default:
;                   Value in conffiles/mass_deriv.conf]
;       BROKFRAC -- FOR BROKEN PL ONLY:  The fraction of the upper
;                   power-law region (in log space) that should be
;                   added to the lower power-law region to ameliorate
;                   small-N issues.  [Default: 0.]
;       SAMPLETD -- Override of sampletd value in mass_deriv.conf, may
;                   be 0 or 1.  [Default: value in mass_deriv.conf]
;
; KEYWORD PARAMETERS:
;       ML      -- Compute the Maximum-Likelihood mass function only.
;                  If set, NREAL is ignored.
;       MEGA    -- Compute a MEGA mass function, with NREAL number of
;                  stacked realizations.  Exclusive of ML -- If both
;                  set, then ML will be chosen.
;       PL      -- Compute a power-law fit to the mass function
;                  [DEFAULT]
;       BROKEN  -- Fit a broken power-law.
;       LN      -- Compute a log-normal fit to the mass function
;       MFNSAVE -- Only compute the desired mass function and save it
;                  to disk WITHOUT fitting a functional form or making
;                  a plot.  This keyword must be used with either ML
;                  or MEGA, but overrides PL and LN.  If ML / MEGA are
;                  not selected, ML will be chosen.
;       VERBOSE -- Cause the routine to contract logorrhoea.
;       SILENT  -- Mum's the word.
;
; OUTPUTS:
;       NONE  [Saves mfn to an IDL save file, and makes plots.]
;
; OPTIONAL OUTPUTS:
;       MFNFN -- Filename of the mass funciton IDL save file
;
; COMMON BLOCKS:
;       OMNI_CONFIG -- The set of configuration structures, read in
;                      from the config files in conffiles/
;       MASS_BLOCK  -- Set of variables read in and computed in
;                      OMNI_MASS_SETUP which are used in other
;                      mass-related routines.
;
; MODIFICATION HISTORY:
;
;       Created:  03/21/13, TPEB -- Initial version -- not complete!!!
;       Modified: 04/02/13, TPEB -- Completed the maximum-likelihood
;                                   section and added the iterative
;                                   version.
;       Modified: 04/03/13, TPEB -- Made memory- and disk-happy by
;                                   removing the randomly-sampled
;                                   distances for each object from the
;                                   ITER structure that is saved to
;                                   disk.  This leaves only the output
;                                   quantities for each iteration in
;                                   the structure.  Also, we now only
;                                   load up 10K distances at a time,
;                                   greatly reducing the amount of
;                                   memory required for the task.
;                                   With this change, very large
;                                   numbers of realizations are
;                                   possible.
;       Modified: 04/03/13, TPEB -- Added functionality to sample dust
;                                   temperature from a Gaussian
;                                   distribution (defined in
;                                   mass_deriv.conf) for each object,
;                                   adding an additional source of
;                                   uncertainty.
;       Modified: 04/15/13, TPEB -- Added functionality to fit a
;                                   broken power-law.
;       Modified: 05/02/13, TPEB -- Added Log-Normal fitting
;                                   functionality, and cleaned up the
;                                   way code sections are parsed and
;                                   run.  Added subroutines for PL and
;                                   LN fitting.
;       Modified: 05/09/13, TPEB -- Brought log-normal code up to
;                                   functional status, added /PL
;                                   keyword to explicitely demand
;                                   power-law fit.
;       Modified: 05/20/13, TPEB -- Added option for naming multiple
;                                   iterations of a given realization
;                                   size in the corresponding IDL save
;                                   filename.
;       Modified: 05/22/13, TPEB -- Added BROKFRAC optional input to
;                                   specify at runtime the fraction of
;                                   overlap the lower PL region has on
;                                   the upper.
;       Modified: 05/29/13, TPEB -- Ensured all output files are
;                                   properly named based on inputs.
;                                   Added SAMPLETD optional input to
;                                   override mconf.sampletd so that
;                                   messing with that config file
;                                   during testing is no longer
;                                   necessary.
;       Modified: 05/30/13, TPEB -- Fixed BROKFRAC bug (typo).
;       Modified: 06/04/13, TPEB -- Added check for PVEC structure IDL
;                                   save file, run DISTANCE_OMNIBUS if
;                                   file is not extant.
;       Modified: 06/10/13, TPEB -- Added MFNSAVE keyword, which
;                                   allows the mass function to be
;                                   saved without going through the
;                                   CPU-gymnastics of computing a
;                                   functional fit to the mfn.
;       Modified: 06/11/13, TPEB -- Changed mfn MEGA structure from
;                                   dist.xx to objs.xx for clarity.
;                                   Added /VERBOSE and /SILENT
;                                   keywords.  Fixed potential bug
;                                   with the use of REFORM. 
;       Modified: 07/17/13, TPEB -- Enforce mass floor as the 90%
;                                   completeness level as computed by
;                                   OMNI_MASS_COMPLETENESS.pro.
;       Modified: 08/21/13, TPEB -- Added support for the new GRS 13CO
;                                   kinematic distance likelihood DPDF
;                                   whereby either KDIST or GRSMATCH
;                                   must be specified for an object to
;                                   be considered for a mass.
;       Modified: 09/03/13, TPEB -- To shrink output plots for MEGA
;                                   cases, set switch to plot line
;                                   rather than open circles for data
;                                   sets larger than N=5000.
;       Modified: 02/04/14, TPEB -- Moved setup section to
;                                   OMNI_MASS_SETUP because of the
;                                   need for duplication in other mass
;                                   routines.  Consequently added the
;                                   MASS_BLOCK COMMON block for ease
;                                   of passing information back from
;                                   the new setup routine.
;       Modified: 02/28/14, TPEB -- Changed loop indices from i -->
;                                   ii, j --> jj for ease of finding.
;                                   Also, modify Monte Carlo mass
;                                   computation to also sample from
;                                   the FLUX DENSITY distribution
;                                   implied by flux +- eflux in the
;                                   catalogs.
;       Modified: 03/13/14, TPEB -- Moved the subroutines OMNI_MFN_PL
;                                   and OMNI_MFN_LN to their own files
;                                   for ease of calling from other
;                                   places.
;
;-

PRO OMNI_MASS_FUNCTION, CONFFILE=cfile, NREAL=nreal, MMAX=m_max,$
                        ML=ml, MEGA=mega, LN=ln, PL=pl, BROKEN=broken, $
                        ITERN=itern, BROKFRAC=brokfrac, SAMPLETD=sampletd, $
                        MFNSAVE=mfnsave, MFNFN=mfnfn, VERBOSE=verbose, $
                        SILENT=silent
  
  COMPILE_OPT IDL2, LOGICAL_PREDICATE
  
  COMMON OMNI_CONFIG, conf, mw, local, dpdfs, ancil, fmt, conffile
  COMMON MASS_BLOCK, pvec, s, iflux, ieflux, m0, hnuk, kappa, nkp
  
  ;; Timing for length of run
  t_start = systime(1)
  
  ;; Parse keywords
  verbose  = KEYWORD_SET(verbose) && ~KEYWORD_SET(silent)
  silent   = KEYWORD_SET(silent)  && ~KEYWORD_SET(verbose)
  nreal    = ~n_elements(nreal) ? 1000L : long(nreal)
  ln       = KEYWORD_SET(ln)        ; Compute a log-normal fit to the mfn?
  pl       = ~ln || KEYWORD_SET(pl) ; Compute a power-law fit to the mfn? [ELSE]
  broken   = KEYWORD_SET(broken)    ; Fit a broken POWER-LAW?
  ml       = KEYWORD_SET(ml)        ; Compute the Maximum-Likelihood mfn?
  mega     = KEYWORD_SET(mega) && ~ml ; Compute a MEGA mfn?
  multi    = ~ml && ~mega      ; Compute multiple realizations of the mfn / fit
  casefit  = pl + 2*ln         ; Add other methods as bitflags (4,8,16,etc)
  itsuf    = n_elements(itern) ? string(itern,format="('_n',I0)") : ''
  brokfrac = n_elements(brokfrac) ? float(brokfrac) : 0.
  mfnsave  = KEYWORD_SET(mfnsave)
  IF mfnsave && multi THEN BEGIN
     ml = 1b
     multi = 0b
  ENDIF
  
  ;; Read in survey-info, galactic-params, & dpdf-params config files
  conf     = omni_load_conf(cfile)
  mw       = omni_read_conffile('./conffiles/galactic_params.conf')
  dpdfs    = omni_read_conffile('./conffiles/dpdf_params.conf')
  mconf    = omni_read_conffile('./conffiles/mass_deriv.conf')
  d        = dindgen(dpdfs.nbins)*dpdfs.binsize + dpdfs.binstart
  sampletd = n_elements(sampletd) ? sampletd GE 1 : mconf.sampletd
  
  status = OMNI_MASS_SETUP(CONFFILE=cfile, VERBOSE=verbose, SILENT=silent, $
                          TYPSTR='mass-function')
  IF ~status THEN RETURN
  
  ;; Read in the SURVEY catalog 
  restore,'./local/'+conf.survey+'_map_locations.sav',VERBOSE=verbose
  
  ;; Limit SURVEY structure to the items saved in the above line
  surv_ind = WHERE_ARRAY(s.cnum, survey.cnum, n_survey)
  IF n_survey NE n_elements(s) THEN message,'Rum now, Think later!!!'
  survey = survey[surv_ind]
  
  
  ;; Parse whether there is a stated maximum mass to be considered
  ;;   KEYWORD passed mmax trumps mconf.mmax
  ;;      If neither is set, mmax = 0., and the maximum is taken as
  ;;      the maximum value of the data (/GETMMAX)
  mmax = KEYWORD_SET(m_max) ? double(m_max) : mconf.mmax
  IF mmax && ~mfnsave THEN BEGIN
     IF ~silent THEN BEGIN
        message,'The maximum mass to be used for fitting functional forms is '+$
                string(mmax,format="(I0)")+' M_sun.',/inf
        print,''
     ENDIF 
     getmmax = 0b
  ENDIF ELSE getmmax = 1b
  
  ;; Use OMNI_MASS_COMPLETENESS to obtain the 90% completeness level
  ;;   for this sample, and use this as the floor for fitting a
  ;;   functional form to the mass function.
  IF ~mfnsave THEN $
     cfn = OMNI_MASS_COMPLETENESS(mscale, ML=ml, MEGA=mega, NREAL=nreal, $
                                  SAMPLETD=sampletd, DMFBINW=binw, /SILENT, $
                                  CLIM90=clim90,CLIM50=clim50,CLIM10=clim10)
  
  
  
  ;;=========================================================================
  ;;=========================================================================
  ;; Choose which kind of mass function to compute
  
  CASE multi OF
     
     ;; Single-shot mass function, whether ML or MEGA.  Includes plotting
     0: BEGIN
        
        ;; Flow chart: ML or MEGA --> provides mfn
        ;; Choose PL or LN --> provides fit information (including plot strings)
        ;; Plot: frame plot w/ CCDF and binned histogram
        ;; Choose PL or LN --> plots fit and information related thereto.
        
        CASE ml + 2*mega OF
           1: BEGIN
              ;; Create the Maximum-Likelihood Mass Function
              
              ndloop = 1
              objs = replicate( {cnum:pvec.cnum,$
                                 dist:dblarr(nkp),$
                                 mass:dblarr(nkp)}, ndloop)
              
              ;; Distance from the PVEC structure
              objs.dist = pvec.stat.duse[0]
              
              ;; Compute the mass function & populate structure
              objs.mass = m0 * s.(iflux) * conf.fluxcor * $
                          (objs.dist / 1.d3) * (objs.dist / 1.d3)
              
              ;; Get entries which conform to the distance region chosen
              dind = where(objs.dist/1.d3 GE mconf.dmin $
                           AND objs.dist/1.d3 LE mconf.dmax, nd)
              darr = (objs.dist)[dind] / 1.d3
              mfn  = (objs.mass)[dind]
              IF ~silent THEN BEGIN
                 ;; print,'DIST: ',m4_stat(objs.dist)/1.d3
                 ;; plothist,objs.dist/1.d3,bin=0.2
                 print,'DIST: ',m4_stat(darr)
                 help,darr
                 print,'MFN: ',m4_stat(mfn)
                 help,mfn
                 
                 message,'For mass-function computation, using the '+$
                         string(nd,format="(I0)")+' objects in the range '+$
                         string(mconf.dmin,format="(F0.1)")+' kpc < d < '+$
                         string(mconf.dmax,format="(F0.1)")+' kpc.',/inf
                 print,'ND: ',nd
              ENDIF
              
              ;; Compute the mass function
              mtot = total(mfn)
              
              ;; Save the mass function to disk
              mfnfn = './masses/save_files/mfn_ml.sav'
              IF verbose THEN message,'Saving to file: '+mfnfn,/inf
              save,objs,mfn,s,iflux,m0,darr,filename=mfnfn,VERBOSE=verbose
              
              ;; ML Mass Function Strings
              plotfn = 'ml_'
              plottit = 'Maximum-Likelihood Mass Function'
              ndfmt = "(I0)"
              
           END
           2: BEGIN
              ;; Create a MEGA Mass Function
              
              ndloop = nreal
              ;; Each element in objs represents ONE LOOP, not ONE OBJECT
              objs = replicate( {cnum:pvec.cnum,$
                                 dist:dblarr(nkp), $
                                 mass:dblarr(nkp), $
                                 flux:dblarr(nkp), $
                                 temp:dblarr(nkp) + mw.td}, ndloop )
              
              IF ~silent THEN $
                 message,'Loading up '+string(ndloop,format="(I0)")+$
                         ' random distances for each object...',/inf
              
              
              ;; Compute temperature distribution, if desired
              IF sampletd THEN BEGIN
                 ;; If sampletd is set, set up what is needed for providing
                 ;;   a random temperature distribution.
                 ;;   Temperature distribution is LOG-NORMAL!!!
                 ;; 0K - 100K base array in 0.01K steps
                 tdarr = dindgen(10001)/100. 
                 lnsig = ( alog(mconf.tdmu + 0.5*mconf.tdfwhm) - $
                           alog(mconf.tdmu - 0.5*mconf.tdfwhm) ) / 2.355
                 tdpdf = gauss_1(alog(tdarr),[1.,alog(mconf.tdmu),lnsig])
              ENDIF
              
              ;; Loop through the keeper objects and populate each with a
              ;;   random sampling from their DPDFs, along with a
              ;;   sample from the FLUX DENSITY array, and
              ;;   TEMPERATURE, if desired.
              FOR ii=0,nkp-1 DO BEGIN
                 
                 ;; Compute the flux density PDF, then load up fluxes!
                 sran = s[ii].(iflux) + [-1.d,1.d]*s[ii].(ieflux)*10 ; 10-sigma
                 dels = sran[1] - sran[0]                    ; Width of above
                 sarr = dindgen(1001)/1000. * dels + sran[0] ; 1000-element arr
                 spdf = gauss_1(sarr,[1.,s[ii].(iflux),s[ii].(ieflux)]) ; PDF
                 sind = where(sarr LT survey[ii].noise, nsi)
                 IF nsi NE 0 THEN spdf[sind] = 0. ; Set prob = 0 for S < noise
                 
                 ;; Load OBJS structure
                 objs.dist[ii] = omni_sample_dpdf(pvec[ii].post,d,ndloop)
                 objs.flux[ii] = omni_sample_dpdf(spdf,sarr,ndloop)
                 IF sampletd THEN $
                    objs.temp[ii] = omni_sample_dpdf(tdpdf,tdarr,ndloop)
              ENDFOR
              
              
              ;; If sampling a dust temperature, generate tempcorr;
              ;;   else set to scalar 1.
              tempcorr = sampletd ? (exp(hnuk/objs.temp)-1.d) / $
                         (exp(hnuk/mw.td)-1.d) : 1.d
              
              ;;***********************************
              ;; Compute the mass function en masse
              objs.mass = m0 * objs.flux * conf.fluxcor * tempcorr * $
                          objs.dist * objs.dist / 1.d6
              
              ;; Get entries which conform to the distance region chosen
              dind = where(objs.dist/1.d3 GE mconf.dmin $
                           AND objs.dist/1.d3 LE mconf.dmax, nd)
              darr = (objs.dist)[dind] / 1.d3
              mfn  = (objs.mass)[dind]
              IF ~silent THEN BEGIN
                 print,'DIST: ',m4_stat(darr)
                 help,darr
                 print,'MFN: ',m4_stat(mfn)
                 help,mfn
              ENDIF
              
              mtot = total(mfn) / float(ndloop)
              nd   = float(nd) /  float(ndloop)
              ndfmt = "(F0.1)"
              
              ;; Save the mass function to disk
              megafn = './masses/save_files/mfn_mega_td'+$
                       string(sampletd,format="(I0)")+'.sav'
              mfnfn = megafn
              IF verbose THEN message,'Saving to file: '+mfnfn,/inf
              save,objs,mfn,s,iflux,m0,darr,filename=megafn,VERBOSE=verbose
              IF ~silent THEN $
                 message,'Computing functional fit for the MEGA mess.  '+$
                         'This may take a while...',/inf
              
              ;; MEGA Mass Function Strings
              plotfn = 'mega_'
              plottit = 'MEGA Mass Function  N = '+string(ndloop,format="(I0)")
              
           END
           ELSE: message,"Don't know what kind of MFN to compute!"
        ENDCASE
        
        ;; If we are only interested in the MFN save file, then RETURN
        IF mfnsave THEN RETURN
        
        ;;----------------------------------------------
        ;; Begin the fitting of a functional form to mfn
        ;; Get type of fit for filename
        CASE casefit OF
           1: BEGIN
              plotfn += 'pl_'
              IF broken THEN $
                 plotfn += 'bf'+string(brokfrac*100.,format="(I02)")+'_'
           END
           2: plotfn += 'ln_'
           3: message,'Error, rummy!'
        ENDCASE
        
        ;;==========================================================
        ;; Set up the generic plot environment, and plot the CCDF
        myps,'./masses/plots/'+plotfn+'mass_fn_td'+$
             string(sampletd,format="(I0)")+'.eps',xsize=6.5,ysize=3.5
        sunsym = cgSymbol('sun')
        
        ccdf = ccdf(mfn)
        Plot,mfn,ccdf,/xlog,xtit='Mass  [M'+sunsym+']',ytit='P(>M)',$
               /nodata,/ylog,xtickformat='exponent10',charsize=1.0,yst=8,$
               ytickformat='exponent10',tit=plottit,xst=5,xmargin=[8,8],$
               xr=[1.,1.d5]
        CASE (n_elements(mfn) GT 5e3) OF ; Make plots faster to render
           0: cgOplot,mfn,ccdf,psym='opencircle',color='blk6',symsize=0.7
           1: cgOplot,mfn,ccdf,color='blk6',thick=15
        ENDCASE
        
        ;;==========================================================
        ;; Now that we have a mass function in the variable mfn, fit
        ;;   something to it, and overplot the fit.
        CASE casefit OF
           
           ;; Fit the power-law slope(s), overplot
           1: BEGIN
              
              alpha = OMNI_MFN_PL(mfn, mmax, mmin, alpha2, mmin2, mmax2, $
                                  BROKEN=broken, BROKFRAC=brokfrac,$
                                  GETMMAX=getmmax, FIXEDMIN=clim90)
              
              ;; Plot the Best-Fit Power-Law Slope & M_min
              qind = where(mfn GE mmin AND mfn LE mmax)
              q    = mfn[qind]
              cf   = (q/mmin)^(1.d - alpha) * ccdf[qind[0]]
              cgOplot,q,cf,thick=5,color='red'
              vline,mmin,/log,color='forest green',thick=3
              vline,mmax,/log,color='cyan',thick=3
              
              ;; Broken power-law
              IF broken THEN BEGIN
                 ;; Plot the Best-Fit Power-Law Slope & M_min
                 qind = where(mfn GE mmin2 AND mfn LE mmax2)
                 q2   = mfn[qind]
                 cf   = (q2/mmin2)^(1.d - alpha2) * ccdf[qind[0]]
                 cgOplot,q2,cf,thick=5,color='orange'
                 vline,mmin2,/log,color='tg2',thick=3
              ENDIF
              
              ;; Draw Legends
              al_legend,/top,/right,/clear,textcolor='forest green',$
                        charsize=0.9,$
                        [cgSymbol('alpha')+' = '+string(alpha,format="(F0.2)"),$
                         'M!dmin!n = '+string(mmin,format="(E0.2)")+$
                         ' M'+sunsym,'N!d>min!n = '+$
                         string(n_elements(q)/float(ndloop),format=ndfmt)]
              IF broken THEN $
                 al_legend,/top,/left,/clear,textcolor='forest green',$
                           charsize=0.9,$
                           [cgSymbol('alpha')+' = '+$
                            string(alpha2,format="(F0.2)"),$
                            'M!dmin!n = '+string(mmin2,format="(E0.2)")+$
                            ' M'+sunsym,'N!d>min!n = '+$
                            string(n_elements(q2)/float(ndloop),format=ndfmt)]
           END
           
           ;; Fit log-normal function, overplot
           2: BEGIN
              
              mu = OMNI_MFN_LN(mfn, sig, mmin, mmax, $
                               GETMMAX=getmmax, FIXEDMIN=clim10)
              
              FORWARD_FUNCTION LNFIT_X
              ;; Plot the Best-Fit Log-Normal and mass limits
              qind = where(mfn GE mmin AND mfn LE mmax)
              q    = mfn[qind]
              cf   = (erfc(LNFIT_X(q,mu,sig))   -erfc(LNFIT_X(mmax,mu,sig)))/$
                     (erfc(LNFIT_X(mmin,mu,sig))-erfc(LNFIT_X(mmax,mu,sig)))*$
                     ccdf[qind[0]]
              cgOplot,q,cf,thick=5,color='red'
              vline,mmin,/log,color='forest green',thick=3
              vline,mmax,/log,color='cyan',thick=3
              
              ;; Draw Legends
              al_legend,/top,/right,/clear,textcolor='forest green',$
                        charsize=0.9,$
                        ['Log-Normal Fit',$
                         cgSymbol('mu')+' = '+string(mu,format="(F0.2)"),$
                         cgSymbol('sigma')+' = '+string(sig,format="(F0.2)")]
              al_legend,/top,/left,/clear,textcolor='forest green',$
                        charsize=0.9,$
                        ['M!dmin!n = '+string(mmin,format="(E0.2)")+$
                         ' M'+sunsym,$
                         'M!dmax!n = '+string(mmax,format="(E0.2)")+$
                         ' M'+sunsym,'N!d>min,<max!n = '+$
                         string(n_elements(q)/float(ndloop),format=ndfmt)]
              
           END
           
           ELSE: message,'Either zero or more than one fit types are '+$
                         'specified.  Taking a shot of rum...'
        ENDCASE
        
        
        ;; Plot the Differential Mass Function
        binw = mconf.dmfbinw NE 0. ? mconf.dmfbinw : fdr_width(alog10(mfn))
        plothist,alog10(mfn),xarr,yarr,bin=binw,/noplot
        yarr = [0,yarr,0] / float(ndloop)
        dx   = mean(xarr[1:*]-xarr[0,*])
        xarr = [xarr[0]-dx,xarr,xarr[n_elements(xarr)-1]+dx]
        
        ;; Use the completeness correction on yarr
        yarr_cor = interpol(cfn, alog10(mscale), xarr)
        ;;yarr /= yarr_cor
        
        ;; Limit to finite values for yarr
        finind   = where(finite(yarr))
        yarr     = yarr[finind]
        xarr     = xarr[finind]
        yarr_cor = yarr_cor[finind]
        
        dmfcol = 'cyan'
        yr = [0.1,10.^(ceil(alog10(max(yarr,/nan))))*10.]
        cgAxis,yaxis=1,yr=yr,/yst,/save,ytit='N per M bin',$
               charsize=1.0,ylog=1,color=dmfcol,ytickformat='exponent10'
        cgOplot,10.^(xarr),yarr,psym=10,color=dmfcol,thick=3
        ;; cgOplot,10.^(xarr),yarr*yarr_cor,psym=10,color='blu4',thick=3
        vline,/log,clim90,thick=5,color='pur4',linestyle=3
        
        
        ;; For the MEGA case, compute and overplot the error bars
        IF mega THEN BEGIN
           IF ~silent THEN $
              message,'Computing differential mass function error bars...',/inf
           nbin = n_elements(xarr)
           hist = fltarr(nreal,nbin)
           ;; For each realization, make distance cuts & histogram
           FOR ii=0, nreal-1 DO BEGIN
              d = objs[ii].dist / 1.d3
              dind = where(d GE mconf.dmin AND d LE mconf.dmax, nd)
              hist[ii,*] = histogram(alog10(objs[ii].mass[dind]), binsize=binw,$
                                     min=min(xarr)-dx/2., nbins=nbin, loc=xh)
           ENDFOR
           ;; Loop through bins to get mu, sig
           sighist = fltarr(nbin)
           FOR jj=0,nbin-1 DO sighist[jj] = stddev(hist[*,jj])
           sighist /= yarr_cor
           sighist = sighist[finind]
           oploterror,10.^xarr,yarr,sighist,errcolor='BLU7',psym=3
           save,objs,mfn,s,iflux,m0,darr,xarr,yarr,sighist,$
                filename=megafn,VERBOSE=verbose
        ENDIF
        
        ;; Plot fit atop the differential mass function
        CASE casefit OF
           ;; Power-law slope(s), overplot
           1: BEGIN
              df = (q/mmin)^(1.d -alpha) * interpol(yarr,10.^(xarr),mmin)
              cgOplot,q,df,thick=5,color='tg3'
              IF broken THEN BEGIN
                 df = (q2/mmin2)^(1.d -alpha2) * interpol(yarr,10.^(xarr),mmin2)
                 cgOplot,q2,df,thick=5,color='forest green'
              ENDIF
           END
           ;; Log-Normal Function
           2: BEGIN
              yf = mpfitpeak(xarr,yarr*yarr_cor,A,nterms=3) ; Get approx norm
              df = gauss_1(alog(q),[A[0],mu,sig])
              ln_cor = interpol(cfn,mscale,q)
              ;; cgOplot,q,df/ln_cor,thick=5,color='TG3'
              cgOplot,q,df,thick=5,color='TG3'
           END
        ENDCASE
        
        ;; Overall stats -- make sure these are computed correctly
        ;;                  for the MEGA variety
        tdstr = sampletd ? 'Log-Normal' : $
                string(mw.td,format = "(I0,' K')")
        al_legend,/bottom,/left,/clear,textcolor='forest green',charsize=0.9,$
                  ['N = '+string(nd,format=ndfmt),$
                   'M!dtot!n = '+string(mtot,format="(E0.2)")+' M'+sunsym,$
                   'T!dd!n = '+tdstr]
        
        ;; Plot x-Axes after the rest to get a clean appearance
        cgAxis,xaxis=0,/xst,/xlog,xtit='Mass  [M'+sunsym+']',$
               xtickformat='exponent10',charsize=1.0
        cgAxis,xaxis=1,/xst,/xlog,xtickformat='blank_axis'
        
        myps,/done
        nreal = 1               ; For output statistics
        
     END
     
     ;;==============================================================
     ;; Do a proper mass function using the full DPDFs
     1: BEGIN
        
        ;;==============================================================
        ;; Do a proper mass function using the full DPDFs
        
        ;; Create structure to hold outputs, based on fit type
        CASE casefit OF
           1: BEGIN             ; Power-Law
              fnsuf = 'pl_'
              str = {alpha:0.d,$
                     mmin:0.,$
                     mmax:0.,$
                     mtot:0.d,$
                     ndist:0L,$
                     nfit:0L}
              IF broken THEN $
                 str = CREATE_STRUCT(str, 'alpha2',0.d, 'mmin2',0., 'mmax2',0.)
           END
           2: BEGIN             ; Log-Normal
              fnsuf = 'ln_'
              str = {mu:0.d,$
                     sig:0.d,$
                     mmin:0.d,$
                     mmax:0.d,$
                     mtot:0.d,$
                     ndist:0L,$
                     nfit:0L}
           END
           3: message,'Error, rummy!  If you got here, '+$
                      'I was drunker than I thought...'
        ENDCASE
        iter = replicate(str, nreal)
        
        ;; Generate the print format code for the status messages
        pfmt = "(I"+string(ceil(alog10(nreal+1)))+")"
        
        ;; Set up the distances structure and loop index
        ndloop = long(nreal < 1e4)
        dloop = -1LL
        objs = replicate({dist:dblarr(nkp)}, ndloop)
        
        ;; If sampletd is set, set up what is needed for providing
        ;;   a random temperature distribution.
        ;;   Temperature distribution is LOG-NORMAL!!!
        IF sampletd THEN BEGIN
           temp = replicate({objs:dblarr(nkp)}, ndloop)
           tdarr = dindgen(10001)/100. ; 0K - 100K base array in 0.01K steps
           lnsig = ( alog(mconf.tdmu + 0.5*mconf.tdfwhm) - $
                     alog(mconf.tdmu - 0.5*mconf.tdfwhm) ) / 2.355
           tdpdf = gauss_1(alog(tdarr),[1.,alog(mconf.tdmu),lnsig])
           ;; cgPlot,tdarr,tdpdf,symcolor='red'
        ENDIF
        
        ;; Now, start looping through the iterations, and make like above
        FOR jj=0LL,nreal-1 DO BEGIN
           
           IF (jj+1) MOD 100 EQ 0 && ~silent THEN BEGIN
              message,'Working realization '+string(jj+1,format=pfmt)+' of '+$
                      string(nreal,format=pfmt)+'...',/inf
              IF nreal GE 1e4 THEN BEGIN
                 etd = t_start + $
                       (systime(1)-t_start) / (double(jj+1) / double(nreal))
                 message,'[ ETD: '+systime(0,etd)+' ]',/inf
              ENDIF
           ENDIF
           
           ;;===================================================================
           ;; Load up 10K random-sample distances at a time.  This
           ;;   DRASTICALLY saves on memory compared to having it in the
           ;;   iter structure, especially since we don't care
           ;;   about the actual distances used for each object for each
           ;;   iteration.
           IF jj MOD ndloop EQ 0 THEN BEGIN
              IF ~silent THEN $
                 message,'Loading up '+string(ndloop,FORMAT="(I0)")+$
                         ' random distances for each object...',/inf
              ;; Loop through the keeper objects and populate each with a
              ;;   random sampling from their DPDFs.
              FOR ii=0,nkp-1 DO $
                 objs.dist[ii] = omni_sample_dpdf(pvec[ii].post,d,ndloop) ;,$
              ;; pvec[ii].stat.kdar,$
              ;; pvec[ii].stat.dtan)
              ;; Also load up random temperatures, if desired
              IF sampletd THEN $
                 FOR ii=0,nkp-1 DO $
                    temp.objs[ii] = omni_sample_dpdf(tdpdf,tdarr,ndloop)
              dloop++           ; Increment the counter
           ENDIF
           ;;===================================================================
           
           ;; Extract distances for this iteration
           darr = objs[jj-(dloop*ndloop)].dist / 1.d3
           dind = where(darr GE mconf.dmin AND darr LE mconf.dmax, nd)
           
           ;; If sampling a dust temperature, generate tempcorr
           IF sampletd THEN BEGIN
              tempcorr = (exp(hnuk/temp[jj-(dloop*ndloop)].objs)-1.d) / $
                         (exp(hnuk/mw.td)-1.d)
              tempcorr = tempcorr[dind] ; use only the DIND objects!!!
           ENDIF ELSE tempcorr = 1.d
           
           
           ;; Compute the mass function & populate ITER structure
           mfn = m0 * s[dind].(iflux) * conf.fluxcor * $
                 darr[dind] * darr[dind] * tempcorr
           iter[jj].mtot = total(mfn)
           iter[jj].ndist = nd
           
           
           ;;===========================================================
           ;; Pick yer fittin' poision
           CASE casefit OF
              1: BEGIN          ; Fit the power-law slope(s)
                 iter[jj].alpha = OMNI_MFN_PL(mfn, mmax, mmin, $
                                             alpha2, mmin2, mmax2, /SILENT, $
                                             BROKEN=broken, BROKFRAC=brokfrac, $
                                             GETMMAX=getmmax, FIXEDMIN=clim90)
                 iter[jj].mmin  = mmin
                 iter[jj].mmax  = mmax
                 q             = mfn[where(mfn GE mmin, nq)]
                 iter[jj].nfit  = nq
                 IF broken THEN BEGIN
                    iter[jj].alpha2 = alpha2
                    iter[jj].mmin2  = mmin2
                    iter[jj].mmax2  = xmax
                 ENDIF
              END
              
              2: BEGIN          ; Fit a log-normal function
                 iter[jj].mu   = OMNI_MFN_LN(mfn, sig, mmin, mmax, $
                                            GETMMAX=getmmax, FIXEDMIN=clim90, $
                                            /SILENT)
                 iter[jj].sig  = sig
                 iter[jj].mmin = mmin
                 iter[jj].mmax = mmax
                 q            = mfn[where(mfn GE mmin AND mfn LE mmax,nq)]
                 iter[jj].nfit = nq
              END
              3: message,'Some other derogatory comment about the programmer.'
           ENDCASE
           
        ENDFOR                  ; End of iteration loop
                
        ;; Save for posterity -- CHANGE THIS TO ./local/ AT SOME POINT!!!
        save,iter,filename='./masses/save_files/iter_mfn_'+fnsuf+'td'+$
             string(sampletd,format="(I0)")+'_'+$
             string(nreal / float(10.^(floor(alog10(nreal)))),$
                    floor(alog10(nreal)),format="(I0,'e',I0)")+itsuf+'.sav',$
             VERBOSE=verbose
        
     END
     ELSE: message,'If this is ever called... Oy... Everclear?',/cont
  ENDCASE
  
  
  ;; Print to screen time statistics
  IF ~silent THEN BEGIN
     time = systime(1) - t_start
     hour = time/3600.
     min  = (time - floor(hour)*3600)/60.
     sec  = time - floor(hour)*3600 - floor(min)*60
     message,'This run had sampletd = '+$
             string(sampletd,format="(I0)")+'.',/inf
     message,'For '+string(nreal,format="(I0)")+' realizations of the mass '+$
             'function, the routine took '+string(floor(hour),format="(I0)")+$
             ' hours, '+string(floor(min),format="(I2)")+' minutes, and '+$
             string(round(sec),format="(I2)")+' seconds.',/inf
  ENDIF
  
  RETURN
END
